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ABSTRACT 

Integral field unit (IFU) data of the iconic Pillars of Creation in M16 are presented. 
The ionisation structure of the pillars was studied in great detail over almost the 
entire visible wavelength range, and maps of the relevant physical parameters, e.g. 
extinction, electron density, electron temperature, line-of-sight velocity of the ionised 
and neutral gas are shown. In agreement with previous authors, we find that the pillar 
tips are being ionised and photo-evaporated by the massive members of the nearby 
cluster NGC 6611. They display a stratified ionisation structure where the emission 
lines peak in a descending order according to their ionisation energies. The IFU data 
allowed us to analyse the kinematics of the photo-evaporative flow in terms of the 
stratified ionisation structure, and we find that, in agreement with simulations, the 
photo-evaporative flow is traced by a blueshift in the position-velocity profile. The 
gas kinematics and ionisation structure have allowed us to produce a sketch of the 
3D geometry of the Pillars, positioning the pillars with respect to the ionising cluster 
stars. We use a novel method to detect a previously unknown bipolar outflow at the 
tip of the middle pillar and suggest that it has an embedded protostar as its driving 
source. Furthermore we identify a candidate outflow in the leftmost pillar. With the 
derived physical parameters and ionic abundances, we estimate a mass loss rate due 
to the photo-evaporative flow of 70 Mq Myr“^ which yields an expected lifetime of 
approximately 3 Myr. 
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1 INTRODUCTION 

Throughout their entire lifetime, massive stars influence 
their immediate surroundings via strong stellar winds, po¬ 
tent ionising radiation and powerful supernovae events. Sim¬ 
ulations show that their feedback is responsible for clearing 
vast bubble-shaped gas voids, inflating expanding HII re¬ 
gions, exposing pillar-like dust and gas lanes and regulating 
the formation of new stellar populations in their surrounding 


natal molecular clouds ([Gritschneder et al.|20I0 

[Peters et al. 

2011 

Bisbas et al.||20II[ [Tremblin et al. [20I2[ 

Walch et al. 

2013 

Dale et al.[[20I3a| [Dale et al.[[20I3b| [Ngoumou et al. 

2015 

1 . The simplified and idealised conditions of the Simula- 


tions can yield important insight on feedback mechanisms, 
but even the state of the art simulations fail to represent the 
complexity of reality, and the details of this feedback have 
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yet to be fully understood. Furthermore, massive star for¬ 
mation dominates the energetics and feedback in star form¬ 
ing galaxies, and properly accounting for the star formation 
feedback is a critical ingredient of galaxy evolution mod¬ 
els. This is the first publication from our FuSIOn (Feedback 
in massive star forming regions: from Simulations to Ob¬ 
servations) study in which we seek to test the predictions 
of the above mentioned numerical simulations, to put bet¬ 
ter constraints on their initial conditions and to ultimately 
understand feedback mechanisms by comparing them with 
observations of high mass stellar feedback. 

Pillar-like structures are a common morphological phe¬ 
nomenon in star forming regions where massive O- and B- 
type stars shape the surrounding material via ionisation and 
stellar winds, and they arise naturally in simulations of clus¬ 
ter formation and evolution that include radiative feedback 


(Gritschneder et al.||20I0 Dale et al.||20I2 Tremblin et al. 
2QI2||Walch et al.||2QI3r Such structures are found to host 
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young stellar objects and are observed to be between the 
hot and hostile environment created by massive stars and 
the material of the parent molecular cloud that has not yet 
been reached by the stellar feedback. The nearby massive 
stars compress, ionise and photo-evaporate the top layers of 
gas that face directly toward them ( Hester|1991 ), and there¬ 
fore ultimately cause the pillars to be destroyed over time. 
Whether indeed the massive stars are able to trigger star 
formation at the pillar-ambient medium interface, and how 
efficient their photo-evaporating effect really is are still open 
questions that need to be addressed. 

One of the prime examples of pillars in a massive star 
forming region is the Pillars of Creation in the Eagle Neb¬ 
ula (M16), which can be seen as protruding from a molecular 
cloud into a vast HII region inffated by the 2-3 Myr old and 
~ 2 kpc distant (Hillenbrand et al. 1993) massive cluster 


NGC 6611, whose brightest and most massive object (the 
03-05 V star W205) lies - in projection - just short of 2 
pc North-East of the Pillars. The Pillars can be divided into 
three main structures, all of them pointing towards the clus¬ 
ter stars of NGC 6611. In this work we will refer to them 
as PI (the easternmost and biggest pillar, subdivided into 
Pla and Plb), P2 (middle pillar) and P3 (the westernmost 
pillar), as is shown in Eig.[^ This pillar-like shape seems to 
have its origin in the fact that the low-density material in 
the bodies of the pillars could be easily removed, but it is 
capped by much denser cores which shield the pillars from 
the intense ionising sources in NGC 6611 ( 01iveira|2008 ). 

These Pillars have been studied by many authors in 
the past two decades since they were first captured in spec¬ 
tacular detail with the Hubble Space Telescope (HST) by 
Hester et al. ( 1996| (hereafter referred to as H96), and as 
it is too great a task to collect and cite all publications re¬ 
garding these famous structures, we can here only give a 
brief historical research overview. H96 first described the 
morphology in great detail and analysed the stratified ion¬ 
isation structure in terms of a photo-evaporative flow and 
identified star forming cometary globules currently emerg¬ 
ing from the Pillars with the HST images obtained in the 
Ha, [SH]A6717,6731 and [OHIJA5007 filters. By comparing 
the optical observations with near-IR images, these authors 
suggest the presence of a population of YSOs in the pillars 
and so called evaporating gaseous globules (EGGs), indi¬ 
cating active star formation in the region. [McCaughrean 
Andersen (2002) however find that only about 15 % of the 
EGGs show signs of active star formation, while the rest 
appear not to host young stellar objects. Eor completeness 
we show the spectrum of an EGG not covered by the HST 
data of H96 which is located to the left of PI (the spectrum 
is shown in Eig. & and its location is shown in Eig. |12^). 
A molecular line and continuum emission study, combined 
with mid-IR observations, was presented by [White et al.| 
(1999), who reported that the pillars are capped by dense 
and relatively cold (~ 10-60 M©) cores in a pre-protostellar 
evolutionary stage. These authors detected no embedded IR 
sources or molecular outflows at the pillar tips, found very 
small (~ 1.7 km/s) velocity gradients along the pillar bod¬ 
ies and estimated the total mass of the Pillars to be ~ 200 
M©. Eukuda et al. ( [2002 ) used millimetre data to analyse 
the age sequence of GO cores, finding a linear age gradi¬ 
ent with more evolved objects closer to NGG 6611 and less 
evolved ones closer to the Pillars, suggesting the propaga¬ 


tion of star formation. Whether this corresponds to triggered 
star formation is debated, as suggested by jlndebetouw et al.j 
(2007), but the presence of young stellar objects (YSOs) and 
water masers near or in the pillar tips is discussed by sev¬ 
eral different authors (Thompson et al. 2002| Healy et al. 


2004 Linsky et al. 2007). In the more recent years, much 


effort has been invested in the simulation of such objects 
to better understand the extent of the stellar feedback, the 
physical effects occurring at the pillar-ambient interface and 


their kinematic structure ( 

Miao et al. 2006 Mizuta et al. 

2006 Ercolano et al.||2012|) 

. To celebrate Hubble’s 25th an- 


niversary, the NASA/ESA/Hubble Heritage Team recently 
released a new, revisited HST image of the Pillars with a 
much higher spatial resolution and a bigger held of view 
than H96, bringing new media attention to these structures. 

As part of the science verification of the optical inte¬ 
gral held spectrograph MUSE at the Very Large Telescope 
(VLT), we obtained a 3x3 aremin mosaic of the Pillars. In 
this work we analyse the ionisation structure of these objects 
in greater detail than ever before, and tie their evapora¬ 
tion and active star formation directly back to the members 
of NGG 6611. Eollowing the theme of our EuSIOn project, 
we seek to compare our observations with the results ob¬ 
tained by [Ercolano et al.|[2012| who created synthetic ob¬ 
servations of pillars under the influence of ionising O- and 
B-type stars in their immediate neighbourhood. They com¬ 
puted line profiles as a function of position within the pillars 
for Ha, [OHI]A5007, [NH]A6584 and [SH]A6717, as well as 
line ratios for the just mentioned emission lines. 

This paper is organised as follows. We introduce both 
the observations and simulations in Section 2; we compute 
the main physical parameters such as electron density and 
temperature, emission line and ratio maps, as well as analyse 
the ionisation structure in Section 3. In Section 4 we describe 
the velocity structure obtained via emission line fitting and 
compute the expected lifetime of the Pillars. Gonclusions are 
presented in Section 5. The appendix is available as online 
supplementary material. 


2 OBSERVATIONS AND SIMULATIONS 
2.1 IFU data 

The optical lEU observations of the Pillars of Greation in 
M16 were obtained with the Multi Unit Spectroscopy Ex¬ 


plorer MUSE instrument (Bacon et al. 2010) mounted on 


the VLT. The data was obtained during the instrument’s 
science verification run in service mode on June 21 and June 
22 (2014) for the proposal 60.A-9309(A) (PI Me Leod), us¬ 
ing the nominal wavelength range spanning from 4750 A to 
9350 A (and a wavelength-dependent spectral resolution of 
~ 75 and 150 km respectively). The Wide Eield Mode of 
MUSE has a held of view of 1x1 aremin^, and we obtained 
two 3x3 aremin^ mosaics with an exposure time of 30 and 
130 seconds each, thus requiring 9 pointings per mosaic. The 
nine fields of the two mosaics are shown in Eig. [^overlaid 
on the HST Ha+[NH] image, and their central coordinates 
and number of spaxels (spectral pixels) are listed in Table 
Each held was observed three times via a dither pattern 
by rotating the instrument position angle by 90 degrees at 
each exposure. The full dataset thus consists of two subsets 
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Table 1. Central coordinates of the MUSE mosaic pointings 
shown in Fig. ^ The number of spaxels (spectral pixels) corre¬ 
sponds to the number of valid spaxels in the Ho; slice of the data 
cube. 


Field no. 

RA (J2000) 

Dec (J2000) 

No. of spaxels 

1 

18:18:51.918 

-13:48:45.90 

94586 

2 

18:18:48.892 

-13:49:23.52 

93936 

3 

18:18:45.866 

-13:50:01.18 

94129 

4 

18:18:48.493 

-13:50:45.96 

94746 

5 

18:18:51.519 

-13:50:08.34 

94723 

6 

18:18:54.545 

-13:49:30.58 

94034 

7 

18:18:57.172 

-13:50:15.50 

94207 

8 

18:18:54.146 

-13:50:53.12 

94690 

9 

18:18:51.120 

-13:51:30.78 

94610 


of 27 exposures each, where the subsets differ in exposure 
time only. The effective spatial resolution of the presented 
data is about 0.2", which corresponds to ~ 1.9x10“^ pc (or 
- 390 AU). 

The 54 exposures were reduced using the standard 
recipes of the MUSE pipeline (vO.18.2) and the required cal¬ 
ibrations. To make the integrated line maps of the entire mo¬ 
saic discussed in Section 3.1, the integrated line maps of the 
single pointings were combined using Montage. After this 
step of making the mosaic a fault in the world coordinate 
system (WCS) of the single fields, resulting from the rotat¬ 
ing dither pattern, manifested itself. This lead to the single 
fields being shifted and rotated with respect to one another, 
and an updating of the WCS of each field by comparing ob¬ 
served stellar coordinates with a catalog was required. This 
was done with the iraf tasks MSCTPEAK and MSCSETWCS 
of the MSCRED package. 

As a counterpart to the iconic and very well known 
HST 3-color composite by [Hester et al. (1996) and the re¬ 
cently released 25th HST anniversary image, the RGB com¬ 
posite (where red is [SH]A6717, green is Ha, and blue is 
[OHIJA5007) of the 30 second exposure is shown in Fig.[^. 
This RGB composite clearly shows the structure and evap¬ 
oration of the Pillars and their tips which are being illu¬ 
minated and ionised by the nearby cluster NGC 6611, and 
also clearly to be seen are the dust lanes at the lower end of 
the Pillars that point towards their natal molecular cloud, 
as well as several detached globules to the left of PI (better 
appreciable in the digital version if this article). 

Figure shows a continuum 3-color composite of the 
region: the pillar tips (PI and P2) are to be seen as reflection 
nebulae illuminated by the nearby stars. The star at the tip 


of P2 corresponds to a known T Tauri star (Sugitani et al. 
2002 Indebetouw et al.||2007 ). 




(b) 


Figure 2. RGB composite (panel a): red is [SIIJA6717, green is 
Ho, and blue is [OIIIJA5007. Continuum RGB composite (panel b, 
red 8050 - 8250 A, green 6850 - 7000 A, blue 5200 - 5400 A). The 
tips of PI and P2 are to be seen as reflection nebulae illuminated 
by the nearby stars. The white circle in panel b indicates the 
continuum source discussed in Section 4.1. 


2.2 Hydrodynamical simulations 

We compare our observations with three-dimensional 
radiation-hydrodynamical simulations performed with the 
smoothed-particle hydrodynamics (SPH) code divine. The 
code is able to model the plane-parallel irradiation by ion¬ 
ising photons of an arbitrary gas distribution, as described 


by Gritschneder et al. (2009), and also includes the effects 


of the diffuse ionising field resulting from recombinations 


(Ercolano Gritschneder|2011 ). 


Ercolano & Gritschneder (2011) simulated the pho- 
toionisation of a (4pc)^ turbulent box initially containing 
474M0 of neutral gas at lOK and with a turbulent RMS 
Mach number of 5. The box is illuminated with a Lyman 
continuum flux of 5x10^ photons cm“^ s~^ and allowed 
to evolve for 500 kyr, corresponding to ~ 0.17 freefall 
times. The turbulent density field is sculpted by the photo- 
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Figure 1. The nine MUSE pointings overlaid on the Hubble Space Telescope Hq:+[NII] image (credit: STScI). Each pointing spans 
approximately 1x1 aremin^. The red lines represent the cuts along which the line intensity profiles were computed (Section 3.2). 



evaporation flows and the pressure in the HII region into a 
complex morphology which includes three prominent pillar 
structures pointing roughly in the direction from which the 
photons are propagating. It is these structures which form 
the basis for our comparative work. 

The SPH simulations and subsequent mocassin (de¬ 
scribed in Section 2.3) calculations do not specifically repro¬ 
duce the pillars in Ml6, as they have smaller mass and the 
radiation field is not as strong as the one produced by NGC 
6611 (see Section 3.3). The comparison we seek to do is more 
of a qualitative one which helps us to see whether the char¬ 
acteristic ionisation structure of the pillar tips is recovered 
in the simulations after running them through the radiative 
transfer calculation, ultimately testing our numerical com¬ 
putations. 


2.3 Monte Carlo radiation-transport calculations 

Since the SPH simulations described above are only 
intended to reproduce the dynamical effects of photoevapo¬ 


ration and ionised gas pressure, they do not lend themselves 
to detailed comparison with observations. For this purpose, 
we post-processed the SPH simulation output with the 


3D Monte Carlo photoionisation code mocassin (Ercolano 
et al. 2003 2005| 2008). mocassin uses a Monte Carlo 


formalism, enabling the selfNconsistent treatment of the 
direct and diffuse radiation field. 

We mapped a 1x3x1 pc cutout of the hydrodynamical 
simulation snapshot onto a uniform 128x384x128 pixel 
grid. We assumed typical HII region gas-phase abundances 
as follows: He/H = 0.1, C/R = 2.2x10“^, N/H = 4.0xl0-^ 
0/H = 3.3x10“^, Ne/H = 5.0xl0“^ S/H = 9.0x10“^ 
The gas is heated mainly by photoionisation of hydrogen 
and other abundant elements. Cas cooling is mainly via 
collisionally excited lines of heavy elements, although con¬ 
tributions from recombination lines, free-bound, free-free 
and two-photon continuum emission are also included. To 
compute the emission line flux from each cell, mocassin 
solves the statistical equilibrium problem for all atoms and 
ions, given the local temperature and ionisation state. 

These mocassin calculations were already used by 
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Ercolano et al. (2012) to create synthetic observations 


of the simulations of Ercolano & Gritschneder (2011). 


Here, for the first time, we make use of the very rich new 
dataset provided by MUSE to compare the hydrodynamic 
simulations in detail to the Pillars of Creation. 


3 IONISATION STRUCTURE 
3.1 Extinction correction 

The extinction toward M16 has been studied by many au¬ 
thors in the past and has been found to vary significantly 


over the members of NGC 6611. Chini & Wargau (1990) 
found that for A < 5500 A the extinction law deviates 


from the standard galactic extinction law (Rv = 3.1, Seaton 
1979), this effect is caused by the dust grains being larger 


than those found in the standard interstellar medium (Or- 
satti et al. 2006). Indeed, Hillenbrand et al. (1993) find a 


typical value of Rv ~ 3.75 for NGC 6611. Therefore, to cor¬ 
rect for extinction along the line of sight toward M16 we 
computed the reddening coefficient C via the Balmer decre¬ 
ment method: 


C = 3.1 X log 


F{Ha) 


- log 


I{Ha)\ 

Km) 


ass uming a reddening curve parametrised by [Cardelli] 


et al. ( 1989| with Rv = 3.1 for A > 5500 A, and Rv 


= 3.75 for shorter wavelengths. An intrinsic flux ratio of 
= 2.86 (corresponding to a case B Balmer 
recombination decrement, [Osterbrock Ferland||2006 ) was 
adopted, and the measured F{Ha)/F{Hl3) flux ratio used 
to compute C was taken to be the mean value obtained 
from the map of that particular line ratio, thus obtaining 
{F{Ha)/F{Hj3)) = 5.93 and therefore C = 0.98. We find 
extinction values between ~ 1.5 and ~ 3.2 mag, which agree 
with literature values ( [Guarcello et al.||2010| . Fig. shows 
co-added spectra of the HII region (top panel), the pillar- 
ambient interface (the pillar tips, middle panel) and an EGG 
(bottom panel). The prominent features in the 8000 A to 
9000 A range are OH sky lines that have not been removed 
from the spectra). Tablelists the detected emission lines 
together with their observed wavelength (column 2), the 
adopted reddening curve (column 4) and the flux ratio prior 
and after correcting for extinction (columns 3 and 5) for the 
example spectrum of the HII region shown in Fig.[^ The in¬ 
tegrated line intensity maps were de-reddened according to 
the appropriate reddening curve value for each wavelength. 
The [OH]A7320,30 as well as the [NHJA5755 have intensities 
(relative to H/3) < 0.07 and were not included in the fit used 
to derive the parameters in Table 

A clear gradient can be seen in the extinction map 
shown in Fig. & where Pla shows a higher overall extinc¬ 
tion, P3 the lowest, and P2 values in between. This is quan¬ 
tified in Fig. Et- where the normalised histograms of circular 
regions (r ~ 3 arcsec, white circles in panel a) in the three 
pillar tips are shown together with the best fit values to each 
histogram. We therefore believe that P3 is the pillar closest 
along the line of sight, behind that lies P2, and behind P2 
lies PI. A description of the 3D geometry of the Pillars is 
given in Section 4. 


3.2 Integrated line maps 

The extinction-corrected integrated line maps were created 
using the moment map subroutine of the python package 
spectral_cubeQ All line and line ratio maps shown in 
this paper (comprised the appendix) are linearly scaled to 
minimum/maximum, please refer to the online version for 
Appendix figures. The Ha, [NHJA6548, [SH]A6717,31 and 
[OHIJA5007 maps are shown in Fig.[^ the maps of the other 
detected lines are shown in Appendix A, available as supple¬ 
mentary online material. H96 discussed the normal-oriented 
striations coming from the pillar-ambient interface in the 
context of photo-evaporation: as the strong ionising radia¬ 
tion from NGC 6611 hits the dense molecular material of the 
pillars, the matter is ionised and the pressure at the pillar- 
ambient matter interface is increased causing matter to flow 
away from the latter, creating a photo-ionised flow. The fact 
that the evaporation is normal to the pillar-ambient inter¬ 
face is an indication that neither the interaction between the 
flow and the ambient matter, nor the interaction of the im¬ 
pinging photons with the outflowing material are able to in¬ 
fluence the direction of the photo evaporative flow (see H96). 
The photo-evaporative striations can be clearly seen in the 
Ha maps (and slightly in the [SHIJA9068 and [ArHI]A7135 
maps). The emission from the other species is either very 
diffuse ([OHI], H/3) or very localised (rest). 

The strongest emission is the Ha emission coming from 
the tip of PI. In general we identify two types of emission: 
a more or less localised line emission that traces the tips of 
the pillars as well as the protrusions on the inner side of PI 
can be seen in the Ha, [Arlll], Hel, [Nil], [01], [OH], [SH], 
[SHI] maps, and a more diffuse emission that surrounds the 
pillars in a halo-like manner ([OHI] and H/3). In the first 
type the emission peak of some lines is sharper and more 
localised than others, e.g. [SH] compared to Ha, and we fur¬ 
thermore see a spatial shift in the peak emission of different 
lines corresponding to a stratified ionisation structure which 
correlates with the ionisation energy of the identified species. 
This stratification was already discussed in H96, who mod¬ 
elled the scenario with CLOUDY84 (see Per land et al.|20l3 
for the latest release of the code) and found a very good 
agreement between the model and the data. While these 
authors were limited to three emission lines (Ha, [SH] and 
[OHI]), MUSE now offers the possibility of investigating this 
stratified ionisation structure in great detail across almost 
the entire optical regime and thus covering almost all of the 
optical ionisation-tracing emission lines. Fig. shows the 
intensity profiles of the emitting line species listed in Table 
l^for the three pillars along the cuts marked in Fig. [^as a 
function of distance from the top of the image (wherever two 
lines of the same ionised state and atom were present these 
were averaged): the stratified ionisation structure is clearly 
seen from the fact that the lines from the ions/atoms with 
higher ionisation energy peak first, followed by the ones with 
lower and lower ionisation energy. The locations of the peaks 
with respect to the [OHI] peak are listed in Table it is 
clear that MUSE is not able to spatially resolve the ionisa¬ 
tion front, as for some lines (especially in P2 and P3 where 
the emission is weaker compared to PI) it is impossible to 
spatially separate the maxima. 


spectral-cube.readthedocs.org 
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Figure 3. Example spectra of the HII region (t op p anel), the pillar-ambient matter interface (middle panel) and the EGG (bottom 
panel) extracted from the positions marked in Fig.|l2^ (see text Section 3.1). The flux is measured in 10“^® ergs s“^ cm“^ pixel” 
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Figure 4. Left: extinction map in terms of (linear scale to maximum/minimum). To help the reader identify the location of the 
Pillars has been marked with a schematic black contour. The white circles in panel (a) correspond to the regions used to quantify the 
extinction for the three Pillars (panel b and text Section 3.1). Right: normalised extinction histograms of circular regions (r ~ 3 arcsec) 
at the tip of the three pillars. 


Table 2. Detected emission lines (column 1) and line centroids obtained via gaussian line fitting (column 2), line fluxes relative to the 
observed flux prior (column 3) and after (column 5) correcting for extinction, and the reddening curve for each line (column 4). 


Line 

Ao6s 

F(A)/F(H/3) 

f(A) 

F(A)/F(H/3)corr 

4861 

4861.42 

1.0T0.019 

0.0 

1.000 

4959 [OIII] 

4959.01 

0.433T0.015 

-0.025 

0.419 

5007 [OIII] 

5006.93 

1.36T0.023 

-0.036 

1.297 

5577 [OI] 

5577.42 

0.298T0.014 

-0.153 

0.244 

5876 Hel 

5875.76 

0.233T0.015 

-0.203 

0.179 

6300 [OI] 

6300.39 

0.382T0.015 

-0.267 

0.269 

6363 [OI] 

6363.86 

0.121T0.014 

-0.276 

0.084 

6548 [Nil] 

6548.12 

0.312T0.015 

-0.303 

0.21 

6563 Ha 

6562.89 

6.476T0.088 

-0.305 

4.343 

6584 [Nil] 

6583.5 

0.975T0.019 

-0.308 

0.552 

6678 Hel 

6678.24 

0.083T0.015 

-0.322 

0.055 

6717 [SII] 

6716.56 

0.197T0.015 

-0.327 

0.129 

6731 [SII] 

6730.94 

0.144T0.015 

-0.329 

0.093 

7065 Hel 

7065.37 

0.047T0.015 

-0.377 

0.029 

7135 [ArIH] 

7135.89 

0.231T0.015 

-0.387 

0.139 

7751 [Arlll] 

7751.19 

0.129T0.014 

-0.475 

0.069 

9068 [SHI] 

9068.69 

0.521T0.015 

-0.531 

0.272 


Fig.j^shows the line intensity profiles obtained from our 
simulations, and although we cannot resolve the separation 
between the Ha, [SII] and [Nil] peaks, the shape and trend 
of the profiles is in good agreement. The distance between 
the [OIII] maximum and the maxima of the other three lines 
is about d 0.015 pc, which agrees with what we find from 
the observations of PI where dHa,[Nii] ~ 0.019 and d[ 5 //] ~ 


0.022. Our models therefore quantitatively predict the ioni¬ 
sation structure of pillar-like objects. 


3.3 Electron temperature and density maps 

The electron density was estimated via the ratio of the 
[SII]A6717 and [SII]A6731 lines. The sulfur ion has a single 
ground level, and two first excited levels are energetically 
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Table 3. Location of the peaks (in 10 ^ pc) of the lines seen in Fig. |^, [^ and[^ with respect to the position of the [Olllj peak (Eio„ 
= 35.12 eV) and their ionisation energy (see text Section 3.2). 



[Arlll] 

Hel 

[Sill] 

Ha,/3 

[Nil] 

[SII] 

[Oil] 

[OI] 

PI 

0.296 

0.739 

1.775 

1.923 

2.218 

2.218 

2.218 

2.514 

P2 

0.296 

0.296 

0.444 

0.444 

0.592 

0.739 

0.739 

1.627 

P3 

0.296 

0.296 

0.296 

0.296 

0.592 

0.592 

0.592 

0.887 

Eion (sV) 

27.63 

24.59 

23.23 

13.59 

14.53 

10.36 

13.62 

- 



RA (J2000) 



RA (J2000) 


(a) 


(b) 



19m00.00s 58.00s 56.00s 54.00s 18hl8m52.00s 


RA (J2000) 



RA (J2000) 


(c) (d) 

Figure 5. Integrated line maps derived from the MUSE data: Ha (a), [NIIJA6548 (b), [SII]A6717,31 (c) and [OIIIJA5007 (d). The flux 
scale bar is in 10“^® erg s“^ cm“^ pixel“^, all maps are are linearly scaled to minimum/maximum. See Section 3.2 for discussion. 
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Figure 6. Line intensity profiles for the emission lines listed in Table for PI, P2 and P3 (panels a, b, and c respectively, see 
text Section 3.2 for discussion). The lines are the same in the first three panels: Ho; (blue), H/3 (yellow), [SIIIJA9068 (solid green), 
[SII]A6717,31 (dashed green), [OIII]A4969,5007 (solid red), [OII]A7320,30 (dashed red), [OI]A5577,6300 (dotted red), Hel A5876,6678 
(magenta), [Arlll]A7135,7751 (cyan), [NII]A6548,6584 (orange). Panel d shows the oxygen lines only for P2: [OIII] (solid line), [Oil] 
(dashed line) and [OI] (dotted line). The profiles in these figures only cover the pillar tips, the do not follow the whole length of the slits 
shown in Fig.^ 


very close and their ratio is thus sensitive to density but 
not temperature. We used the analytical solution given in 


McCall (1984) for a three-level atom: 


[5'//]A6717 

[5'//]A6731 


R[sii] — 1-49 


1 + 3.77x 
1 + 12.Sx 


( 1 ) 


where x = 0.01Ne/T^^‘^ and assuming that T = 10^ K, 
so that the electron density Ne in cm“^ is given by 


Ne = 


R[sii] — 1-49 
5.6713 - 12.8R[sii] 


X 10^ 


( 2 ) 


The tips of the pillars show the highest densities with 
values > 2000 cm“^ and we find decreasing density values 
along the pillar bodies, while the surrounding gas is at den¬ 
sities < 250 cm“^ (Fig.[^). The pillar tips having the high¬ 
est electron densities is a consequence of them also having a 


much higher column density than the pillar bodies (White 
et al.|19^ ): the dense tips act like caps that shield the pillar 
bodies from the stellar feedback. Fig. shows the electron 
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Figure 7. Line intensity profiles for the emission lines obtained from our simulations (panels a and b) along the (horizontal) slits shown 
in the Ho; intensity map (panel c). The colour-coding is the same as in Fig. See text Section 3.2 for discussion. 


density map of our simulations, derived via Eq. 2. The den¬ 
sity of the HII region is comparable to what we observe, as 
is the fact that the exposed pillar protrusions that are being 
ionised display the highest density values. These are lower 
by a factor of ~ 2, which is consistent with the fact that 
the radiation field implemented in the simulations is weaker 
than what is observed in M16. Another factor adding to this 
is that the simulated pillars, while roughly the same size as 
the M16 Pillars, are only about 1/lOth in the mass. Fig. 

shows that the histogram of the electron density follows 
a lognormal distribution with a slight tail toward the low- 
density regime. The lognormal density distribution is well 
recovered in our simulations as is shown in Fig. [^, where 
the electron density as measured from the SPH grid (real) 
is plotted together with the electron density derived via the 
[SII] line ratio (derived). In the real electron density his¬ 
togram it appears that in the simulations the high and low 


density values are both washed out by the line-ratio calcu¬ 
lation, effect arising from smoothing along the line of sight. 
A lognormal electron density distribution has been observed 
for the ionised interstellar medium ( Redheld Falcon|2Q08 


Hill et al.|2007 ) where the driver for the shape is isothermal 


turbulence. In the case of the Pillars of Creation we suggest 
that shocks could produce such a distribution. The presence 
of shocks could in principle be inferred from emission line 
ratios, but, as we show in Section 3.5, these shock indica¬ 
tors may be washed out - and therefore not detected - by 
the dominant inhuence of photoionisation. 


The electron temperature is usually determined via the 
forbidden line ratio ([OIII]A5007+[OIII]A4959)/[OIII]A4363, 
but the [OIII]A4363 is not covered by MUSE and we there¬ 
fore derived the electron temperature Te via the [Nil] dou¬ 
blet according to equation 5.5 in [Osterbrock fc Ferland] 
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(2006) (assuming that the term in the denominator is neg¬ 
ligible) 




2.5 X 10^ 
ln{0.164R) 


(3) 


where R = [iV77]A6548,8584/[iV7/]A5755. From the 
electron temperature map shown in Fig. it is clear that 
the high density regions at the tip of the pillars and the arc¬ 
like structures in PI display lower temperatures with respect 
to the surrounding ambient gas, although the map presents 
a certain level of noise coming from the weak [Nil]5755 line. 


3.4 Abundance tracers and ionic abundances 

Among abundance parameters the most used are the oxy¬ 
gen R23 and the sulfur S23 parameters, which are defined 
as R23 = ([OII1A3726,3729+[OIII1A4959,5007)/H/? (|pSid 
et al.|1979|| and S 23 = ([SI I]A6717,31+[SIII]A9068,9532)/H/3 
( Vilchez fc Esteban||199^ , and have been extensively used 
to determine galactic as well as extragalactic chemical abun¬ 
dances, since these differ as a function of position within a 
galaxy and can therefore be used to study star formation his¬ 
tories and evolutionary scenarios. Other authors have used 
these parameters, together with the excitation measured by 
the ratio of [OII]/[OIII], to study the ionisation structure 
of HII regions ( [Garcia- Benito et al. |2010| [Monreal-Ibero 
et al. dere, we exploit these parameters to distinguish 

between the various regions in our data (pillars, surround¬ 
ing HII region, stars, young stellar objects, and protostellar 
outflows). This is demonstrated in Fig. where S23 (de¬ 
fined as S23 = ([SII]A6717,31+[SIII]A9068)/H/3 because of 
the lacking coverage of the [SIII]A9532 line) is plotted ver¬ 
sus [OII]/[OIII] (panel a) for a sub-region that highlights the 
tip of P2 (panel b): the clearly recognisable different popula¬ 
tions in the scatter plot correspond to ambient matter (red), 
the pillar-ambient interface (green), the pillar matter (blue), 
a known T Tauri star (magenta, see Section 4.1) and a sulfur- 
abundant region (orange) of panel b. The brushed S23 map 
of the entire mosaic is shown in Fig. the stars have been 
intentionally left white. Fitting in the picture where Pla 
lies behind NGC 6611 and the other two pillars along the 
line of sight (see Section 4.2), we see the pillar body of Pla 
(marked in blue) contaminated by the HII region (marked 
in red). The sulfur-abundant region is clearly visible in the 
S23 parameter map (Appendix A) and is visible as a slight 
over-density of the order of ~ 2000 cm“^ in the electron 
density map. As will be discussed in Section 4.1, we suggest 
that this region corresponds to where the blue-shifted lobe 
of a bipolar jet originating from a deeply embedded proto¬ 
star is pushing its way out through the pillar material. The 
fact that it is only visible in the integrated intensity maps 
of the [SH] lines and in [OI]A6300 indicates that the outflow 
is not energetic enough to ionise the detected atoms with 
Eion ^ Eion,[sii] = 10.36 eV. 

Furthermore, we derived ionic and total abundances 
of the HII region and the pillar tips via the I RAF task 
ionic. For this, we produced co-added spectra of a re¬ 
gion representative of the HII region matter and a region 
representative of the pillar tips (see Fig. |^, which were 
fitted with a gaussian routine to determine emission line 


intensities. The electron density was computed via equa¬ 
tion 2, and this was then used to derive ionic abundances 
and the electron temperature with I RAF. The derived ionic 
and total abundances are listed in Table The oxygen 
ionic abundance ratios 0^/H+ and 0^^/H+ were obtained 
with the [OH]A7320,30 and [OHI]A4959,5007 lines respec¬ 
tively and assuming T([OH]) T([NH]) and T([OHI]) 

T([SHI]), as we do not have the wavelength coverage for 
the [OH]A3727,29 and the [OHI]A4363 lines to determine 
T([OHI]) and T([OH]). The total S/H abundance was de¬ 
termined by taking into account the unobservable ionisation 
stages and therefore us ing the ionisation c orrection factor 
(IGF) for as in Hagele et al. (2008). As is expected 

from the fact that molecular material is being ionised at the 
pillar tips, these present the higher ionic abundance values 
when compared with the HII region. 


3.5 Line ratio maps and BPT diagrams 


BPT ( Baldwin et al.|p-981| ) diagrams compare collisionally 
excited lines like [OHIJA5007, [NH]A6584 and [SH]A6717, 31 
to hydrogen recombination lines (Ho, H/3) and use the for¬ 
bidden/recombination line ratios as indicators of the num¬ 
ber of ionisations per unit volume. The position of an ob¬ 
ject on a BPT diagram is used for classification purposes, 
e.g. to distinguish star forming galaxies from active galac¬ 
tic nuclei (AGN), as gas excited by photoionisation rather 
than shock ionisation occupies different regions in such a di¬ 
agram. iKewfe^^eT^T (2001) combined photoionisation and 
stellar population synthesis models to put an upper limit 
on the location of star-forming sources in the BPT diagram. 
This upper limit is given by a maximum starburst line (la¬ 
belled as KeOl in Fig. 13 from Kewley et al.||2001 ) of pure 
stellar photoionisation models, and galaxies lying above this 
line are considered AGN dominated because for them to be 
found here an additional excitation mechanism other than 


stellar photoionisation is required. Kauffmann et al. (2003) 


modified the KeOl plot to account for galaxies that have 
contribution from both star formation and AGN (labelled 
as Ka03 in Fig. 13 from Kauffmann et al.|2003 ). 


The position on the BPT diagram as a diagnostic tool is 
widely used for spatially unresolved observations of galax¬ 
ies or nebulae where the integrated emission of the whole 
source is used to compute the emission line ratio. BPT 
diagrams have been used for IFU observations of HII re¬ 


gions ( Monreal-Ibero et al.|20ir Garcfa-Benito et al.|2010 ) 
as well, but since these are spatially resolved observations 
caution is advised: one is comparing a 3D structure with 
a ID photoionisation model, and in a 2D projection of a 
3D structure the line ratio will show variations from pixel 
to pixel because of the local gas conditions along the line 
of sight on each pixel (jErcolano et al.|2012|). Maps of the 
[NH]A6584/Ha, [SH]AA6717,6731/Ha and [OHI]A5007/H/3 
line ratios are shown in Fig. |12[ 

Here we investigate the commonly used line ratios as 
indicators of the ionisation structure: [OHI]A5007/H/3 is a 
tracer of the highly ionised gas and [NH]A6584 /Hq; as well as 
[SH]AA6717,6731/H(a are tracers of the less ionised gas. The 
[Nil]/Ha and [SH]/Ha line ratio maps are quite similar and 
show significantly higher line ratios at the tip of the pillars as 
well as on the edges of the wave-like structures in the pillars’ 
main bodies, correlating with the position of the peaks of the 
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Figure 8. The electron density map derived via Eq. |^(a), the electron temperature derived via Eq. |^(b), and the electron density of 
the simulation computed with the sulfur line ratio (c) (see text Section 3.3). 
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Table 4. Derived ionic and total abundances of O and S (see Section 3.4). 


Region 

12+log(0+/H+) 

12+log(02+/H+) 

12+log(0/H) 

12+log(S+/H+) 

12+log(S2+/H+) 

12+log(S/H) 

HII 

8.00 

8.16 

8.39 

6.05 

7.53 

7.59 

Pillar tip 

8.76 

8.22 

8.87 

6.79 

7.65 

7.77 


30000 

27000 

24000 

21000 

18000 

15000 

12000 

9000 

6000 

3000 


Electron temperature (K) 
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log(Electron number density (cm-’) logfElectron number density) (cm-) 

(a) (b) 

Figure 9. Histogram of the electron density of the observed data (panel a) and the simulated pillars (panel b). See text Section 3.3 for 
discussion. 


localised emission; the [OIII]/H/3 ratio on the other hand 
has its lowest values at the positions of the [SII]/H(a and 
[Nil]/Ha peaks, as one would expect in an ionisation picture 
where the double ionised elements are found in the less dense 
material more exposed to the ionising radiation rather than 
in the dense molecular environment shielded from the stellar 
influence. 


Fig. shows the [OIII]A5007/Hy5 versus 

[NII]A6584/Ha and [SII]AA6717,6731/Ha diagnostic 
diagrams for the entire observed region (panel a) and for 
a circular region containing the tip of P2 (panel b), the 
lines correspond to the maximum starburst lines from 
Kewley et al. (2001) and Kauffmann et al. ( 2003| ) (to better 


identify structures in panel a, panel c shows a normalised 
2D histogram of the panel a). From these diagrams it 
is clear that the line ratios lie in the the region below 
KeOl and therefore in the region where HII regions are 
expected to be found. Also taking into account the fact 
that if shocks were strongly contributing to the excitation 
one would expect [Nil]/Ha > -0.1 and [SH]/Ha > -0.4 
(Allen et al. 2008[ ) but we do not see values that high, we 
suggest that the dominant excitation mechanism at the 
pillar tips is photoionisation. This agrees very well with the 
BPT diagrams of our SPH simulations (see Fig. 4 and 5 in 
Ercolano et al.[|2012 ), where the majority of data points lie 


in the photoionisation region. The simulated BPT diagrams 
(Figures 4 and 5 in Ercolano et al.|[?012 ) do however show 
a small percentage of points lying in the shock-dominated 
region, but as the simulations do not include shock models 
these line ratios are due to 3D effects. Fig. shows the 
BPT diagram of the tip of P2 color coded according to the 
regions identified in Fig, p^. 


4 EMISSION LINE FITTING 

MUSE’s detectors undersample the line spread function, and 
we therefore proceeded by first normalising and stacking sev¬ 
eral emission lines in each single pixel (thus assuming that all 
the lines in one pixel originate from the same emitting gas) 
and then by fitting the resulting velocity spectrum in each 
pixel with a gaussian function to produce a velocity map[^ 
The lines used for the stacked spectrum are Ha, the two [Nil] 
lines, the two [SH] lines, [OI]A6300,6363 and HeIA6678. The 
undersampling of single lines is shown in Fig. |15[ where the 
stacked spectrum and spectra of single lines (Ha, [NH]A6548 
and [NH]A6584 in panel a, H/3 and [SHI]A9068 in panel b) 
are shown, together with the resolution at 4600 A (Av = 150 
km s“^) and at 9300 A (Av = 75 km s“^). The early ver¬ 
sion of the data reduction pipeline we used did not include 
a radial velocity correction, and the derived velocities are 
therefore topocentric. Furthermore, during the observations 
the sky lines at A5577 A and A6300 A are calibrated with 
a set of arc lamps, calibration that yields an offset between 
the arc and observed lines; this offset is then propagated 
to the rest of the spectrum, but in a linear way such that 
the wavelength dependency of the correction is lost. This 
leads to a velocity offset of the red (e.g. [SHI]) and blue 
(e.g. [OHI], H/3) lines with respect to the lines close to the 
calibration wavelength which can be seen in Fig.[^(a) and 
(b), where the topocentric velocities of single emission line 
species are plotted along the slits shown in Fig. Further¬ 
more, when computing local standard of rest velocities from 


^ The PYTHON packages spectral_cube (spectral- 
cube.readthedocs.org) and PYSPECKiT ( [Ginsburg &: Mirocha| 


20111 were used for the fitting routine. 
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Figure 10. O 23 = [OII]/[OIII] vs. S 23 for the tip of P2 (panel a) and the S 23 parameter map of the same region (panel b). The colours 
in panel a highlight different populations in the scatter plot of panel a which correspond to the same coloured regions in panel (b), as is 
discussed in Section 3.4. 


the MUSE data with I RAF and comparing them to veloc¬ 
ities obtained form CO observations ( PoundpQQS ylsr ~ 
20 - 20 km s“^), we find that the two are offset by about 
20 km s“^, offset we cannot attribute to a physical origin. 
Further investigation of the science verification data and the 
early version of the used data reduction pipeline is needed. 
Because of these problems we only report relative velocities: 
for each line emitting species the mean velocity of the HII 
region was computed, value which was then subtracted to 
the whole velocity map of that line. The relative velocities 
reported are therefore velocities relative to the HII region. 


4.1 Velocity structure 


Based on CO J=3-2 data White et al. (1999) report complex 


velocity fields within the pillars and a systematic large scale 
velocity gradient of about 1.7 km s~^ along P2. The resolu¬ 
tion of MUSE does not allow the detection of such a small 
velocity gradient, but Fig. [T^ shows the velocity map ob¬ 
tained with the above mentioned emission line fitting. The 
pillars are clearly distinguishable from the ambient gas, and 
indeed the pillars are blue-shifted with respect to the sur¬ 
rounding gas because of the photo-evaporative flow along 
the line of sight that is moving radially away from the pil¬ 
lar surfaces. This means that, when moving from the pillar 
bodies facing our way along the line of sight to the pillar- 
ambinet interface, we see the pillar border redshifted with 
respect to the pillar body as a geometrical consequence of 
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Figure 11. S 23 parameter map of the entire mosaic, the colour-coding is the same as in Fig. |10| 



\ / 


Figure 17. Sketch of the normal velocity vector, corresponding 
to the photo-evaporative flow, to demonstrate that the pillar tips 
appear redshifted with respect to the main pillar bodies because 
of the geometrical effect of the vector component along the line 
of sight (see text Section 4.1). 


the normal vector moving out of the line of sight (this is 
seen best in P3, and is illustrated in Fig. 17). The velocities 
of the Pillars and the HII region were estimated first by ex¬ 
tracting circular regions for each of them via brushin^{Fig. 
|16^ ) , and by then fitting the velocity histograms of these re¬ 
gions with Gaussian distributions (Fig. [20|. They are listed 
in Table m 

The sulfur-abundant region identified at the tip of P2 
via the S 23 parameter is blueshifted {vrei ~ -20 km s“^) with 
respect to the main pillar body (v^eZ ~ -7 km s“^), as can 
be seen in Fig, and at the very tip of the same pillar a 
redshifted region is identified {vrei ~ 10 km s“^). The red 


^ The technique of brushing^ also known as graphical exploratory 
data analysis allows the user to manually select speciflc data 
points or subsets from an image or a plot by interactively drawing 
regions on the latter two. 


and blue-shifted regions can be connected with a line with 
a position angle of approximately 107.3 degrees, on which a 
source identified as a candidate T Tauri star with a disk-like 
structure ( Sugitani et al.||2002 Sugitani et al.|[2Q07 


M16 ES-2 in Thompson et al.|2002 ), as well as water maser 
(Healy et al. 2004) are found. We suggest that the blue- 
and redshifted regions trace the two counterparts of a bipo¬ 
lar outflow originating from an accreting protostar or young 
stellar object. As discussed by Healy et al. (2004), the water 
maser is too far away from the candidate T Tauri star to be 
excited by the latter and the two are therefore not associ¬ 
ated. With these observations we are unable to say whether 
the bipolar outflow originates from the T Tauri star or from 
a deeply embedded object, but considering the fact that the 
T Tauri star seems to lie somewhat off the hypothetical di¬ 
rection of the jet, we tentatively suggest that the latter is the 
case. Healy et al. (2004) suggest that the excitation source of 


the maser is a Class 0 protostar due to the lack of a detected 
near-IR source, but high angular resolution sub-mm obser¬ 
vations are required to identify the maser driving protostar. 
Optical protostellar jets generally have observe d radial ve- 
locities from a few tens up to hundreds km s~^ (Bally et al. 


2007), implying that the one detected here is a low-velocity 
outflow with \rei,jet ~ 20 km s“^. The blue-shifted counter¬ 
part is tunnelling its way through the pillar material and is 
only now starting to emerge, and the redshifted counterpart 
on the other hand is flowing out from the pillar tip but is im¬ 
mediately dispersed in the strong ionising and unfavourable 
conditions of the HII region. 


We also identify a candidate outflow on the inner side of 
PI, which is shown in Fig. A10 of the appendix. Also in this 
case we identify both a blueshifted and a redshifted conn- 
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Figure 12. Line ratio maps linearly scaled to maximum/minimum: [Nil]A6584 /Hq; (a), [SII]AA6717,6731 /Hq: (b) and [OIII]A5007/H^ 
(c). The white circles in panel a indicate the regions where the co-added spectra were extracted. The leftmost circle also corresponds to 
the location of the EGG discussed in Section 1. 


terpart in the velocity map that have relative velocities of 
~ -10 km and ~ 16 km respectively. This candidate 
outflow does not appear to correspond to a sulfur-abundant 
region and is not seen in the S23 map. If this is indeed an 
outflow, we suggest that the absence of the sulfur-indicator 
comes from the fact that the blue lobe is only just starting to 
emerge from the pillar body and did not yet have the time 

,[sii] = 10.36 e V. [Healy 


to ionise atoms with E^, 




et al. 


et al. 


(2004) do not detect water masers in PI, but Sugitani 


(2002) identify candidate T Tauri star right outside 


the pillar which corresponds to a continuum source in Fig. 

n 


The presence of ongoing star formation is also observed 
in our SPH simulations (Gritschneder et al. (2010) and Fig. 
A9), where 4 objects have formed at the pillar tips and are 
still present 0.5 Myr into the simulation. These objects all 
have angles of 20-50 degrees between the disk (or accretion 
feature) and the direction of radiation (which in the simula¬ 
tions is roughly perpendicular to the long axis of the pillars). 
|Raga et~ar (2010) simulated outflows from stars forming at 
the tips of pillars and find that the orientation of the out¬ 
flows depends on the orientation of the pillars along the line 
of sight. For pillars that are seen face on, the projected out¬ 
flow orientations are found to be perpendicular to the main 


log([SII]/HQ) 































































The Pillars of Creation revisited with MUSE 17 




log([NII]A6584/HQ) 



-1.6-1.4-1.2-1.0-0.8-0.6-0.4-0.2 0.0 
log([SII]A6717,31/HQ) 


I I I , 

0 10 20 30 40 50 60 70 80 


0 8 16 24 32 40 48 56 64 


(c) 

Figure 13. Diagnostic BPT diagrams for the entire observed area (panel a) and the tip of P2 only (panel b). Panel c shows a normalised 
2D histogram of panel a. See text Section 3.5 for discussion. The KeOl and Ka03 lines separate regions where photoionisation (below the 
lines) and AGN feedback (above the lines) dominate. Ka03 accounts for a mixed feedback of both AGN and photoionisation. 


pillar axes. Then again, Smith et al. (2010) report that in 
their HST Herbig-Haro jet sample the outflows are not pref¬ 
erentially found to be perpendicular to the pillars. Here, the 
projected outflow orientation with respect to the P.A. of ~ 
137° of P2 is ~ 30°, which corresponds to an angle of only ~ 
7° between outflow and the direction of radiation. The disk 
of the jet-driving source is thus supposedly perpendicular 
to the incident radiation, which contradicts what is found 
from the simulations, where the accretion disks of the form¬ 
ing sources are preferentially aligned with the direction of 
radiation (Gritschneder, in preparation). 


Knowing that the ionisation structure of the photo- 
evaporative flow is (spatially) stratified, the question 
whether these emission lines are also offset in position- 
velocity space occurs naturally. Fig. (panels c and d) 
shows the relative velocit^j^of the single atomic species (ob¬ 
tained by simultaneously fitting all the detected lines for 
each species) of P2 and P3 along the same slits used for 
the intensity profiles shown in Fig.[^ The pillar-ambient in¬ 
terface is clearly recognisable from the dip in velocity of ~ 

relative to the HII region 
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Figure 14. BPT diagram of the tip of P2. The colours in panel (a) highlight different populations in the BPT diagram which correspond 
to the same coloured regions in panel (b). 
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Figure 15. Stacked spectrum (red in both panels), shown together with the Ha (green circles), [NIIJA6548 (blue squares) and [NIIJA6584 
(yellow triangles) spectra (panel a), as well as the H/3 (blue squares) and [SIIIJA9068 (yellow triangles) (panel b). See Section 4 for 
discussion. 


7 and 10 km s~^ for P3 and P2 respectively, the dip to¬ 
ward bluer velocities occurring because we are probing the 
radially outward directed photo-evaporative flow which is 
thus moving toward us along the line of sight. The neutral 
[01] line does not follow the same rise-and-fall profile as the 
ionisation-tracing lines, while The [Oil] lines are not shown, 
as their weak intensity to produce a reliable velocity map. 
Panels (a) and (b) of Fig, show the uncorrected (topocen- 
tric) line of sight velocity: the redshift of the [Sill] line and 
the blueshift of the H/3 and [OIII] lines are due to the wave¬ 
length calibration in the MUSE pipeline, as described at the 
beginning of this section. Fig, shows the analogous plot 
for a synthetic pillar: the profile seen in the observations is 


well recovered by our simulations, but the magnitude of the 
blueshift is of the order of about 3-5 km s“^. As already 
discussed, we don’t expect for the simulations and the ob¬ 
servations to agree quantitatively because of the different 
physical conditions. 


4.2 Geometry 

Another interesting open question is the 3D geometry of the 
Pillars and their spatial orientation with respect to NGC 
6611, and because for the first time we have the velocity in¬ 
formation combined with an extinction map and integrated 
line intensities for the entire structures, this can now be ad- 
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(b) 


Figure 16. Map of relative velocity (to the the HII region) of the entire mosaic (panel a) and a zoom-in on P2 (panels 2). The coloured 
circles in panel a correspond to the regions used for the histograms in Fig. |20| The black line in panel b corresponds to a hypothetical 
orientation of the bipolar outflow with position angle 107.3 degrees, the black circle is centred on the coordinates of the candidate T 
Tauri star, while the white circle corresponds to the water maser detected in March 2002 by |Healy et al.|2004| (see text Section 4.1). 
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Figure 18. Velocity profile of P2 (panels a and c) and P3 (panels b and d) along the same slits as Fig. The lines are: Ha (blue), 
H/3 (yellow), [SIII]9068 (solid green), [SII]6717,31 (dashed green), [OIII]4969,5007 (solid red), [OI]5577,6300 (dotted red), Hel 5876,6678 
(magenta), [Arlll]7135,7751 (cyan), [Nll]6548,6584 (orange). Panels (a) and (b) show the observed line of sight velocity (not corrected 
for radial velocity) with the wavelength-dependent offset between the red ([Sill]), blue ([Olll], H/3) and green (rest) lines originating 
from the wavelength calibration applied during the observations (see text). Panels (c) and (d) are the same as the other two, but the 
velocities are now relative to the mean velocity of the region between 0 and 0.2 pc (corresponding to the Hll region). See text Section 
4.1 for discussion. 


dressed. In our attempt to understand the geometry of the 
Pillars, we consider the following: 

• from the integrated intensity and electron density maps 
it is clear that ionisation mainly occurs at the tip of PI and 
P2, that P3 seems to be less exposed to the ionising radiation 
on the side we can see, and that the exposed protrusions 
along PI are also being (though somewhat less than the 
tip) ionised 

• the extinction map (Fig. [^) shows a gradient from PI 
(most extinction) to P3 (least extinction) 

• the tips of PI and P2 are seen as a reflection nebulae in 
the continuum 3-color composite (Fig. [^, while P3 is more 
seen as a silhouette illuminated from behind 


• the three pillars have different velocities, as is shown in 
Fig. [2^ and Table 

• the velocity of the lower part of PI is consistent with 
the velocity of P2 and P3, and seems inconsistent with the 
upper part of PI 


With molecular line data. Pound (1998) find that PI is 
in reality composed of two separate parts, the upper part 
(PIa) being behind the ionising stars of NGC 6611 with a 
tail inclined away from us and the lower part (Plb) being 
in front of them with its tail is inclined toward us. P2 is 
a single structure with its tail pointing toward us, while 
for P3 their data is inconsistent and inconclusive. From our 
observations it is clear that Pla is the least blue-shifted and 
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Figure 19. Velocity profile of a simulated pillar along the same 
slit as Fig. 



(km/s) 

Figure 20. The histograms of the velocity (relative to the HII 
region, see Section 4.1) of the pillars and the HII region obtained 
from circular regions extracted from the velocity map. 


the most extincted of the three pillars, and its tip shows 
the highest degree of ionisation. Together with the positive 
velocity gradient reported in Pound (1998) and the fact that 
the protrusions along its body show weak signs of ionisation, 
we confirm that Pla is behind the ionising O stars and its 
tip is pointing toward us along the line of sight, while Plb 
is in front of the stars and inclined with its tip pointing 
away from us. Because of the lack of ionisation along their 
bodies, we suggest that both P2 and P3 are also in front of 
the ionising sources and that both are inclined with their 
tips pointing away and their tails pointing toward us, that 
P3 is closer to the observer but because it is slightly less 
blue-shifted and its tip shows less ionisation than P2, it is 
inclined with its tip pointing toward us. We sketch this in 
Fig. 1^ to make it easier for the reader to visualise. 


Table 5. Line of sight velocities relative to the HII region derived 
from fitting gaussian distributions to the histograms in Fig. |20| 
(see Section 4.1). 



Pla 

Plb 

P2 

P3 

Vrel (km S ^) 

-3.43 

-5.34 

-7.14 

-6.62 


4.3 Mass loss rate and lifetime 


Because forbidden lines such as [OI]A6300 or [SII]A6731 are 
optically thin they can be used to estimate the mass of the 
line emitting matter, since the line luminosity is proportional 
to the number of line emitting atoms along the line of sight. 
For the pillar-ambient interface region (Te ~ 11870 K, Ne ~ 
2800 cm“^, Fig.[^, the ionisation fraction is Uon ~ 0.01 and 
the critical density is Nc,[s//]A 673 i — 3.9x10^ cm“^, and 
therefore N>Nc. Because the signal-to-noise ratio is much 
higher for the [SII] line than for the [01] line, we use the 
former to derive the mass loss rate of the pillars. This line 
originates from the transition from the ^ 03/2 (level 1 ) to 
the ^ 83/2 (level 2 ) transition which occurs w ith a rate A 21 
= 8.8x10“^ s“^ ( Osterbrock <V Ferland|2006 ). Following the 
method described in |Hartigan et al.| (1995), we compute the 
mass of the line emitting matter with expressions for the 
luminosity Lerai of the [SII]A6731 line and the total number 
of sulfur atoms in the lower level, 771 


Lers, = fm e-p(-^)(l+f) 


m = 


Vi 

r vis) 1 

\V{HP 



Tftot 


Mu 


limH 


which correspond to equations A 2 and A5 in |Hartigan| 
et ah ( 1995| ) and where gi = g 2 = 4 are the statistical weights 
of the two levels, hi>i 2 = 2.95 xl0“^^ erg is the energy of the 
transition, ? 7 (S)/? 7 (H) = 5.88x10“^ (Table [^, 771 / 77 ( 6 ') ~ 
0.9 (computed with the I RAF task ionic), and assuming 
77(H)/77tot ~ 0.921 and ji = 1.24 ( Allen|1973t . With an elec¬ 
tron density of 2800 cm“^ and a temperature of about 10^ K 
characteristic for the pillar tips, the mass of the line emitting 
matter then becomes 


2^69X10- (l + l) 

The mass loss rate can then be computed via M = 
Mv//, where v is the velocity of the photo evaporative out¬ 
flow and I the average size of the line emitting region. In 
order to estimate the total luminosity of the [SII]A6731 line 
we subtracted the stellar continuum in the [SII]A6731 inte¬ 
grated intensity map and summed the line intensity of all 
pixels to obtain Lersi — 130 L©. With v 
the velocity map in Fig. 16 


8 km s ^ from 


and I = 3 xlO^® cm (measured 
from the SII intensity map), M ^ 70 M© Myr“^. With a to¬ 


tal mass of the pillars of Mtot ~ 200 M© ( White et al.|1999t , 
we estimate the lifetime of the pillars to be approximately 
3 Myr. This however is an upper limit, as the [SII] line may 
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Figure 21. Sketch of the 3D geometry of the Pillars of Creation with respect to the ionising stars in NGC 6611 (not to scale). See text 
Section 4.2. 


not be tracing the bulk flow and we therefore might be un¬ 
derestimating the mass loss rate. Pound ( |1998| ) estimated 
the expected lifetime of the Pillars to be about 20 Myr, but 
as suggested by Westmoquette et al.||2013 (who find a life¬ 
time of about 2 Myr for pillars near NGC 3603) we believe 
that the used value of the curvature radius was too small, 
so they overestimated the lifetime. 

We suggest that the identified outflow will not con¬ 
tribute significantly to the mass loss rate. Indeed, we com¬ 
pute the mass of the outflow to be M = 0.29 M© Myr“^ 
(with Ne ~ 1400 cm“^, Lerai — 0.18 L©), which is only ~ 
0.4 % of the total mass loss rate. 


5 CONCLUSIONS 


For the first time the Pillars of Creation were imaged using 
optical integral field spectroscopy. This technique allowed us 
to map the entire structures and obtain medium resolution 
spectra for each pixel in the 3x3 aremin^ mosaic, which were 
fitted with Gaussian profiles to produce integrated line and 
line of sight velocity maps in many spectral lines. Also, we 
computed physical parameters such as the electron density, 
the electron temperature and the extinction for the entire 
mosaic. 


In agreement with Pound (1998), we find that the three 
Pillars are actually 4 distinct structures that lie at different 
distances along the line of sight and have different inclination 
angles, some pointing away and others pointing toward the 
observer. The consequence of their position and inclination 
is a different degree of ionisation not only of the pillar tips, 
but also of the pillar bodies. We therefore see the strongest 
emission from ionised species from the tip of Plb, which is 
pointing toward us but lies behind the ionising O stars of 
NGC 6611 and we therefore see weak emission coming from 
the exposed structures along its body as well. This we do 
not see in Plb, P2 and P3, as these all lie in front of the 


ionising stars, and furthermore for Plb and P3 because they 
are inclined with their tips facing the observer and their tails 
pointing away from us. The effect of their different inclina¬ 
tion angles is that the 4 structures show different velocities: 
all of them are blue-shifted with respect to the ambient HII 
region because we are tracing the photo-evaporative flow, 
but because this flow is normal to their surfaces, we will see 
them at different velocities. 


The kinematic analysis also allowed us to reveal a pos¬ 
sible protostellar outflow at the tip of P2. Of this outflow, 
we can identify both lobes as a blue and a redshifted coun¬ 
terpart. The jet seems to be almost parallel to the incident 
ionising radiation, in fact it is inclined by only 7 degrees 
with respect to the pillar body and the radiation. This does 
not seem to agree with what we find in our smoothed par¬ 
ticle hydrodynamic simulations, in which we also find stars 
forming at the pillar tips, but their accretion disks are pref¬ 
erentially aligned with the direction of radiation. We explain 
this by arguing that the angle between the outflow and the 
pillar bodies is the consequence of projection effects, as is 
discussed in Raga et al. (2010). The outflow is clearly traced 
by the S23 parameter, which, if plotted against an indicator 
of the degree of ionisation such as [Oil]/[GUI], can be used 
in future studies with IFU data of star forming regions to de¬ 
tect outflows in dense molecular material. In our data this is 
the only outflow we can detect in both the velocity map and 
via the S23 parameter. From the kinematics we also identify 
a candidate outflow on the inner side of PI, which however 
does not show a sulfur abundance. Using IFU data it would 
be interesting to apply this method to other star forming 
regions that have an embedded young stellar population to 
look for outflows that have not yet fully emerged from the 
molecular material and would therefore have been missed in 
previous studies. Furthermore, the outflow and its driving 
source be an interesting target for the Atacama Large Mil¬ 
limeter/submillimeter Array, which has the required high 
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sensitivity and the frequency coverage to analyse the veloc¬ 
ity field and the outflow in great detail. 

The Pillars show a classical stratified ionisation struc¬ 
ture in which line emitting species peak in a sequential man¬ 
ner throughout the ionisation interface of the molecular ma¬ 
terial, in agreement with previous authors. Typical electron 
densities at the pillar tips are > 2000 cm“^, whereas the 
electron temperature in this region is in the 8000 K - 10 000 
K range. The comparison with our SPH+radiative trans¬ 
fer computations is encouraging, as we are able to repro¬ 
duce the ionisation structure, the kinematics and physical 
parameters. Finally, by computing the mass of the [SII] line- 
emitting gas, we estimate the mass loss rate to be of about 
70 M© Myr“^, which, assuming the total mass of the Pil¬ 
lars to be ~ 200 M©, yields an expected lifetime of about 3 
Myr, timescale which does not depend of the presence of an 
outflow at the tip of the middle pillar. 
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APPENDIX A: 


This appendix is available as supplementary online material. 
Here we show the extinction corrected integrated line maps 


of the detected emission lines (Fig. |Al| to Fig. A8r, linearly 
scaled to minimum/maximum, the flux is in units of 10~^° 
erg s~^ cm~^ pixeU^), the S 23 parameter as well as a disk¬ 
like structure in one of our simulated pillars (Fig. A9) and 
a zoom-in of the velocity map onto the tip of PI where the 
candidate outflow is located (Fig. AlO). 


This paper has been typeset from a JJ]X/ DTeX file prepared 
by the author. 
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Figure Al. NII6584, SIII9068. 
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Figure A2. 016300, 015577. To help the reader identify the location of the Pillars has been marked with a schematic black contour. 
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(a) 



(b) 

Figure A3. SII6717, SII6731. 
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Figure A4. ArIII7135, ArIII7751. 
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Figure A5. HeI5876, HeI6678. 
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Figure A6. HeI7065, 016363.To help the reader identify the location of the Pillars has been marked with a schematic black contour. 
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Figure A7. OII7320, OII7330. 
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Figure A8. H/5, S 23 parameter (see text for explanation). 
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Figure A9. Zoom-in onto one of the simulated pillars: a disk-like structure is seen tracing a forming stellar object at the tip of the pillar 
(see text Section 4.1). 
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Figure AlO. Zoom-in of the velocity map shown in Fig. 16a onto PI. The black line shows the orientation of the outflow candidate 
discussed in Section 4.1. 
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